function X = simu_brown(X0, v, T, sw)
% Simulate Brownian motion
%
%   X = simu_brown(X0, v, T);
%   X = simu_brown(X0, v, T, sw);
%       simulates Brownian motions.
%
%       Inputs:
%       - X0:       the initial locations, a d x n matrix with each
%                   column corresponding to one point
%       - v:        the variance parameter
%       - T:        the simulation length
%       - sw:       smooth weight.
%
%       Outputs:
%       - X:        the simulated results in d x n x (T+1) matrix, each
%                   page corresponds to a time point. In particular, 
%                   X(:,:,1) equals X0.
%

% Created by Dahua Lin, on Nov 24, 2010
%

%% verify input arguments

if ~(isfloat(X0) && ndims(X0) == 2)
    error('simu_brown:invalidarg', 'X0 should be a numeric matrix.');
end

if ~(isfloat(v) && isscalar(v) && v >= 0)
    error('simu_brown:invalidarg', 'v should be a non-negative real scalar.');
end

if ~(isnumeric(T) && isscalar(T) && T == fix(T) && T >= 0)
    error('simu_brown:invalidarg', 'T should be a non-negative integer scalar.');
end

if nargin < 4
    sw = 0;
else
    if ~(isfloat(sw) && isscalar(sw) && sw >= 0)
        error('simu_brown:invalidarg', 'sw should be a non-negative real scalar.');
    end
end
    


%% main

% special cases

if T == 0
    X = X0;
    return;
end

if v == 0
    X = repmat(X0, [1, 1, T+1]);
    return;
end

% general case

noise = randn([size(X0) T+1]) * v;
X = bsxfun(@plus, X0, cumsum(noise, 3));

% smoothing

[gs,gt,gw] = gridgraph1d(T+1, sw);
g = gr_edgelist.from_edges('u', T+1, gs, gt, gw);

Y = reshape(permute(X, [3 1 2]), T+1, numel(X0));
Y = laplacesm(g, 1, Y);
X = permute(reshape(Y, [T+1 size(X0)]), [2 3 1]);

